import processing.core.*; import processing.opengl.*; import java.applet.*; import java.awt.*; import java.awt.image.*; import java.awt.event.*; import java.io.*; import java.net.*; import java.text.*; import java.util.*; import java.util.zip.*; import javax.sound.midi.*; import javax.sound.midi.spi.*; import javax.sound.sampled.*; import javax.sound.sampled.spi.*; import java.util.regex.*; import javax.xml.parsers.*; import javax.xml.transform.*; import javax.xml.transform.dom.*; import javax.xml.transform.sax.*; import javax.xml.transform.stream.*; import org.xml.sax.*; import org.xml.sax.ext.*; import org.xml.sax.helpers.*; public class electronsphere extends PApplet {

int frame = 0;

Electronsphere es;

public void setup(){
  background(50,50,50);
  size(400,400,P3D);
  //initialize global
  global = new Global();
  global.init();
  //initialize electron sphere
  es = new Electronsphere();
}

public void draw(){
  render();
}

public void render(){
  frame++;
  background(50,50,50);
  pushMatrix();
  translate(width/2,height/2,0);
  stroke(255);
  fill(255,50);
  es.rotate((mouseY-height/2)*.00007f,-(mouseX-width/2)*.00007f);
  es.step();
  es.render();
  popMatrix();
}

public float randomize(float $n,float $r){
  return $n+random(-$r/2,$r/2);
}

class Electron extends Point{
  float fx,fy,fz;
  int fn;
  Vector v;
  boolean locked = false;
  int faces = 0;
  Electron(float $x,float $y,float $z){
    super($x,$y,$z);
    v = new Vector($x,$y,$z);
    reset_forces();
  }
  public void repel(Electron[] $electrons){
    for(int i=0;i<$electrons.length;i++){
      if($electrons[i]!=this){
        float d = get_distance_to($electrons[i].x,$electrons[i].y,$electrons[i].z);
        if(d>0){
          float f = force/sq(d);
          float dx = $electrons[i].x-x;
          float dy = $electrons[i].y-y;
          float dz = $electrons[i].z-z;
          if(!$electrons[i].locked){
            $electrons[i].add_force((-f/2)*(dx/d),(-f/2)*(dy/d),(-f/2)*(dz/d));
          }
          if(!locked){
            add_force((f/2)*(dx/d),(f/2)*(dy/d),(f/2)*(dz/d));
          }
        }
      }
    }
  }
  public void add_force(float $fx,float $fy,float $fz){
    fx += $fx; fy += $fy; fz += $fz; fn++;
  }
  public void average_forces(){
    if(fn==0){ return; }
    fx /= fn; fy /= fn; fz /= fn;
  }
  public void reset_forces(){
    fx = 0; fy = 0; fz = 0; fn = 0;
  }
  public void step(){
    average_forces();
    v.vx = x+fx; v.vy = y+fy; v.vz = z+fz;
    v.unitize();
    v.scale(radius);
    move_to(v.vx,v.vy,v.vz);
    reset_forces();
  }
  public void render(){
    stroke(255,100,100,100);
    v.render(zero_point);
    stroke(255);
    super.render();
  }
}

float radius = 150; //sphere radius
float force = -2000000; //the force with which electrons repel eachother
int tnelectrons = 100; //the desired number of electrons
int tnfaces = tnelectrons*2-4; //the total number of faces
int tiling_start = 900; //the number of frames the electrons have to resolve
int tiling_end = 2700; //give up on tiling after this frame

class Electronsphere{
  Electron[] electrons = new Electron[tnelectrons];
  Face[] faces = new Face[10000];
  int nelectrons = 0;
  int nfaces = 0;
  Electronsphere(){
    for(int i=0;i<tnelectrons;i++){
      electrons[nelectrons] = new Electron(random(-radius,radius),random(-radius,radius),random(-radius,radius));
      nelectrons++;
    }
  }
  public void generate_face(){
    if(nfaces>=tnfaces){ return; }
    Electron[] ns = {electrons[floor(random(nelectrons-.0001f))],electrons[floor(random(nelectrons-.0001f))],electrons[floor(random(nelectrons-.0001f))]};
    Electron[] ons = new Electron[3];
    int tries = 0;
    while(tries<10){
      tries++;
      arraycopy(ns,ons);
      Face f = new Face(ns);
      Point cp = f.get_circumcenter_point();
      float cm = mag(cp.x,cp.y,cp.z);
      if(cm<distance_tolerance){ return; }
      cp.x *= radius/cm; cp.y *= radius/cm; cp.z *= radius/cm;
      ns = get_closest_electrons(cp.x,cp.y,cp.z,3);
      if((ons[0]==ns[0]||ons[0]==ns[1]||ons[0]==ns[2])&&(ons[1]==ns[0]||ons[1]==ns[1]||ons[1]==ns[2])&&(ons[2]==ns[0]||ons[2]==ns[1]||ons[2]==ns[2])){
        f = new Face(ns);
        if(!f.is_acute()){ return; }
        for(int i=0;i<nfaces;i++){
          if((faces[i].p1==ns[0]||faces[i].p1==ns[1]||faces[i].p1==ns[2])&&(faces[i].p2==ns[0]||faces[i].p2==ns[1]||faces[i].p2==ns[2])&&(faces[i].p3==ns[0]||faces[i].p3==ns[1]||faces[i].p3==ns[2])){ return; }
        }
        f.orient(zero_point);
        faces[nfaces] = f;
        nfaces++;
        return;
      }
    }
  }
  public Electron[] get_closest_electrons(float $x,float $y,float $z,int $n){
    Electron[] temp = new Electron[nelectrons];
    float[] ds = new float[nelectrons];
    //get distances to each Electron
    for(int i=0;i<nelectrons;i++){
      ds[i] = dist(electrons[i].x,electrons[i].y,electrons[i].z,$x,$y,$z);
      temp[i] = electrons[i];
    }
    //bubble sort according to distance
    for(int i=0;i<nelectrons-1;i++){
      for(int j=0;j<nelectrons-1-i;j++){
        if(ds[j]>ds[j+1]){
          Electron temp1 = temp[j]; temp[j] = temp[j+1]; temp[j+1] = temp1;
          float ds1 = ds[j]; ds[j] = ds[j+1]; ds[j+1] = ds1;
        }
      }
    }
    //chop it to match number desired
    temp = (Electron[]) subset(temp,0,$n);
    return temp;
  }
  public void rotate(float $rx,float $ry){
    float cosrx = cos($rx); float sinrx = sin($rx);
    float cosry = cos($ry); float sinry = sin($ry);
    for(int i=0;i<nelectrons;i++){
      electrons[i].rotate(cosrx,sinrx,"x",0,0,0);
      electrons[i].rotate(cosry,sinry,"y",0,0,0);
    }
  }
  public void step(){
    if(frame>tiling_start&&frame<tiling_end){
      generate_face();
    }    
  }
  public void render(){
    if(frame<tiling_start){
      for(int i=0;i<nelectrons-1;i++){
        electrons[i].repel((Electron[]) subset(electrons,i+1,nelectrons-i-1));
      }
    }
    for(int i=0;i<nelectrons;i++){
      electrons[i].step();
      electrons[i].render();
    }
    for(int i=0;i<nfaces;i++){
      faces[i].render(true);
    }
  }
}

float distance_tolerance = .0001f;
float slope_tolerance = .001f;
float angle_tolerance = .001f;
float point_speed = .1f;
Point zero_point;
Global global;

//------ GLOBAL ------//
//the global class holds all other default classes
//its only purpose is for organization and to imitate actionscript
class Global{
  int npoints,nsegments,nvectors,nlines,npolygons,npolylines,nplanes,nfaces;
  Point[] points = new Point[1000000]; Segment[] segments = new Segment[1000000]; Vector[] vectors = new Vector[1000000]; Line[] lines = new Line[1000000]; Polygon[] polygons = new Polygon[100000]; Polyline[] polylines = new Polyline[100000]; Plane[] planes = new Plane[100000]; Face[] faces = new Face[100000];
  Polygon stage;
  Global(){
    npoints = nsegments = nvectors = nlines = npolygons = npolylines = nplanes = nfaces = 0;
  }
  public void init(){
    Point[] ps = new Point[4]; ps[0] = new Point(-width/2,-height/2,0); ps[1] = new Point(width/2,-height/2,0); ps[2] = new Point(width/2,height/2,0); ps[3] = new Point(-width/2,height/2,0);
    stage = new Polygon(ps);
    zero_point = new Point(0,0,0);
  }
  public void echo(){
    println("---- GLOBAL ----");
    println("Points: "+npoints);
    println("Segments: "+nsegments);
    println("Vectors: "+nvectors);
    println("Lines: "+nlines);
    println("Polygons: "+npolygons);
    println("Polylines: "+npolylines);
    println("Planes: "+nplanes);
    println("Faces: "+nfaces);
  }
  public void render(int $n){
    if($n>=128){
      for(int i=0;i<nfaces;i++){
        fill(255,240,240);
        stroke(255);
        strokeWeight(1);
        faces[i].render();
      }
      $n -= 128;
    }
    if($n>=64){
      for(int i=0;i<npolylines;i++){
        fill(255,250,250);
        stroke(0);
        strokeWeight(1);
        polylines[i].render();
      }
      $n -= 64;
    }
    if($n>=32){
      for(int i=0;i<npolygons;i++){
        fill(255,250,250);
        stroke(0);
        strokeWeight(1);
        polygons[i].render();
      }
      $n -= 32;
    }
    if($n>=16){
      noFill();
      stroke(0,0,255);
      strokeWeight(1);
      for(int i=0;i<nlines;i++){
        lines[i].render();
      }
      $n -= 16;
    }
    if($n>=8){
      noFill();
      stroke(0,255,255);
      strokeWeight(1);
      for(int i=0;i<nvectors;i++){
        vectors[i].render(zero_point);
      }
      $n -= 8;
    }
    if($n>=4){
      noFill();
      stroke(255,0,255);
      strokeWeight(1);
      for(int i=0;i<nsegments;i++){
        segments[i].render();
      }
      $n -= 4;
    }
    if($n>=2){
      fill(0);
      stroke(0);
      strokeWeight(1);
      for(int i=0;i<npoints;i++){
        points[i].render();
      }
      $n -= 2;
    }
    if($n>=1){
      noFill();
      stroke(255,0,0);
      strokeWeight(1);
      stage.render();
      $n -= 1;
    }
  }
}

//------ POINT ------//
class Point{
  int n;
  float x,y,z,ox,oy,oz,oox,ooy,ooz,tx,ty,tz,otx,oty,otz;
  Point(float $x,float $y,float $z){
    n = global.npoints;
    global.points[global.npoints++] = this;
    x = ox = oox = tx = otx = $x;
    y = oy = ooy = ty = oty = $y;
    z = oz = ooz = tz = otz = $z;
  }
  public void reset(){
    x = oox; y = ooy; z = ooz; ox = oox; oy = ooy; oz = ooz; tx = oox; ty = ooy; tz = ooz; otx = oox; oty = ooy; otz = ooz;
  }
  public void match(Point $p){
    x = $p.x; y = $p.y; z = $p.z; ox = $p.ox; oy = $p.oy; oz = $p.oz; tx = $p.tx; ty = $p.ty; tz = $p.tz; otx = $p.otx;  oty = $p.oty; otz = $p.otz;
  }
  public void move(float $dx,float $dy,float $dz){
    ox = x; oy = y; oz = z; x += $dx; y += $dy; z += $dz; tx = x; ty = y; tz = z;
  }
  public void move_to(float $x,float $y,float $z){
    ox = x; oy = y; oz = z; x = $x; y = $y; z = $z; tx = x; ty = y; tz = z;
  }
  public void direct(float $dx,float $dy,float $dz){
    otx = tx; oty = ty; otz = tz; tx += $dx; ty += $dy; tz += $dz;
  }
  public void direct_to(float $tx,float $ty,float $tz){
    otx = tx; oty = ty; otz = tz; tx = $tx; ty = $ty; tz = $tz;
  }
  public void rotate(float $cosa,float $sina,String $axis,float $x,float $y,float $z){
    float dx,dy,dz;
    if($axis=="x"){
      dz = z-$z; dy = y-$y;
      z = $z+dz*$cosa-dy*$sina; y = $y+dy*$cosa+dz*$sina;
    }else if($axis=="y"){
      dx = x-$x; dz = z-$z;
      x = $x+dx*$cosa-dz*$sina; z = $z+dz*$cosa+dx*$sina;
    }else if($axis=="z"){
      dx = x-$x; dy = y-$y;
      x = $x+dx*$cosa-dy*$sina; y = $y+dy*$cosa+dx*$sina;
    }
    tx = x; ty = y; tz = z; ox = x; oy = y; oz = z; otx = x; oty = y; otz = z;
  }
  public float get_distance_to(float $x,float $y,float $z){
    return dist(x,y,z,$x,$y,$z);
  }
  public void step(){
    ox = x;
    oy = y;
    oz = z;
    x += (tx-x)*point_speed;
    y += (ty-y)*point_speed;
    z += (tz-z)*point_speed;
  }
  public void echo(int $indent){
    String indent = "";
    while(indent.length()<$indent){
      indent += "  ";
    }
    println(indent+"---- POINT #"+n+" ----");
    println(indent+"x: "+x+"  tx: "+tx);
    println(indent+"y: "+y+"  ty: "+ty);
    println(indent+"z: "+z+"  tz: "+tz);
  }
  public void render(){
    beginShape();
    vertex(x-2,y-2,z);
    vertex(x+2,y-2,z);
    vertex(x+2,y+2,z);
    vertex(x-2,y+2,z);
    endShape(CLOSE);
  }
}

  //------ VECTOR ------//
class Vector{
  int n;
  float vx,vy,vz,m;
  Vector(float $vx,float $vy,float $vz){
    n = global.nvectors;
    global.vectors[global.nvectors++] = this;
    vx = $vx;
    vy = $vy;
    vz = $vz;
    m = get_magnitude();
  }
  Vector(Point $p){
    n = global.nvectors;
    global.vectors[global.nvectors++] = this;
    vx = $p.x;
    vy = $p.y;
    vz = $p.z;
    m = get_magnitude();
  }
  public float get_magnitude(){
    return mag(vx,vy,vz);
  }
  public void unitize(){
    m = get_magnitude();
    if(m==0){ return; }
    vx = vx/m; vy = vy/m; vz = vz/m;
  }
  public void add(Vector $v){
    vx += $v.vx; vy += $v.vy; vz += $v.vz;
  }
  public void subtract(Vector $v){
    vx -= $v.vx; vy -= $v.vy; vz -= $v.vz;
  }
  public void scale(float $s){
    vx *= $s; vy *= $s; vz *= $s;
  }
  public void echo(int $indent){
    String indent = "";
    while(indent.length()<$indent){
      indent += "  ";
    }
    println(indent+"---- VECTOR #"+n+" ----");
    println(indent+"vx: "+vx);
    println(indent+"vy: "+vy);
    println(indent+"vz: "+vz);
  }
  public void render(Point $p){
    beginShape();
    vertex($p.x,$p.y,$p.z);
    vertex($p.x+vx,$p.y+vy,$p.z+vz);
    endShape();
  }
}

public float dot_product(Vector $v1,Vector $v2){
  return ($v1.vx*$v2.vx)+($v1.vy*$v2.vy)+($v1.vz*$v2.vz);
}
public float dot_product_1p(Vector $v1,Point $p2){
  return ($v1.vx*$p2.x)+($v1.vy*$p2.y)+($v1.vz*$p2.z);
}
public float dot_product_2p(Point $p1,Point $p2){
  return ($p1.x*$p2.x)+($p1.y*$p2.y)+($p1.z*$p2.z);
}

public Vector cross_product(Vector $v1,Vector $v2){
  return new Vector($v1.vy*$v2.vz-$v1.vz*$v2.vy,$v1.vz*$v2.vx-$v1.vx*$v2.vz,$v1.vx*$v2.vy-$v1.vy*$v2.vx);
}

//------ SEGMENT ------//
class Segment{
  int n;
  float l;
  Point p1,p2;
  Segment(Point $p1,Point $p2){
    p1 = $p1; p2 = $p2;
    n = global.nsegments;
    global.segments[global.nsegments++] = this;
    l = get_length();
  }
  public void reset(){
    p1.reset();
    p2.reset();
  }
  public void match(Segment $s){
    p1.match($s.p1);
    p2.match($s.p2);
  }
  public float get_length(){
    return dist(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z);
  }
  public float get_slope(){
    //xy coordinates only
    if(p2.x!=p1.x){
      return ((p2.y-p1.y)/(p2.x-p1.x));
    }else{
      return 999999999;
    }
  }
  public boolean is_ray_intersecting(float $x,float $y,boolean $inclusive){
    //xy coordinates only :: for inside-outside testing, etc. :: ray points right, starting from the given coordinates
    float pbua1 = (p2.x-p1.x)*($y-p1.y)-(p2.y-p1.y)*($x-p1.x);
    float pbu2 = p2.y-p1.y;
    float pbub1 = $y-p1.y;
    if(pbu2==0){ return false; }
    float pbua = pbua1/pbu2;
    float pbub = pbub1/pbu2;
    if($inclusive){
      if(pbua>=0&&pbub>=0&&pbub<=1){ return true; }
    }else{
      if(pbua>0&&pbub>0&&pbub<1){ return true; }
    }
    return false;
  }
  public boolean is_line_intersecting(Line $l,boolean $inclusive){
    //xy coordinates only
    float pbua1 = (p2.x-p1.x)*($l.s.p1.y-p1.y)-(p2.y-p1.y)*($l.s.p1.x-p1.x);
    float pbu2 = (p2.y-p1.y)*($l.s.p2.x-$l.s.p1.x)-(p2.x-p1.x)*($l.s.p2.y-$l.s.p1.y);
    float pbub1 = ($l.s.p2.x-$l.s.p1.x)*($l.s.p1.y-p1.y)-($l.s.p2.y-$l.s.p1.y)*($l.s.p1.x-p1.x);
    if(pbu2==0){ return false; }
    float pbua = pbua1/pbu2;
    float pbub = pbub1/pbu2;
    if($inclusive){
      if(pbub>=0&&pbub<=1){ return true; }
    }else{
      if(pbub>0&&pbub<1){ return true; }
    }
    return false;
  }
  public boolean is_segment_intersecting(Segment $s,boolean $inclusive){
    //xy coordinates only
    float pbua1 = (p2.x-p1.x)*($s.p1.y-p1.y)-(p2.y-p1.y)*($s.p1.x-p1.x);
    float pbu2 = (p2.y-p1.y)*($s.p2.x-$s.p1.x)-(p2.x-p1.x)*($s.p2.y-$s.p1.y);
    float pbub1 = ($s.p2.x-$s.p1.x)*($s.p1.y-p1.y)-($s.p2.y-$s.p1.y)*($s.p1.x-p1.x);
    if(pbu2==0){ return false; }
    float pbua = pbua1/pbu2;
    float pbub = pbub1/pbu2;
    if($inclusive){
      if(pbub>=0&&pbub<=1&&pbua>=0&&pbua<=1){ return true; }
    }else{
      if(pbub>0&&pbub<1&&pbua>0&&pbua<1){ return true; }
    }
    return false;
  }
  public boolean is_coord_on(float $x,float $y){
    //xy coordinates only
    float d = get_distance_to($x,$y);
    if(d<=distance_tolerance){
      return true;
    }
    return false;
  }
  public float get_distance_to(float $x,float $y){
    //xy coordinates only
    float u = (($x-p1.x)*(p2.x-p1.x)+($y-p1.y)*(p2.y-p1.y))/sq(dist(p1.x,p1.y,p2.x,p2.y));
    float x = p1.x+u*(p2.x-p1.x);
    float y = p1.y+u*(p2.y-p1.y);
    if(x>p1.x&&x>p2.x){ x = max(p1.x,p2.x); }
    if(x<p1.x&&x<p2.x){ x = min(p1.x,p2.x); }
    if(y>p1.y&&y>p2.y){ y = max(p1.y,p2.y); }
    if(y<p1.y&&y<p2.y){ y = min(p1.y,p2.y); }
    return dist(x,y,$x,$y);
  }
  public Point get_closest_point(float $x,float $y){
    //xy coordinates only
    float u = (($x-p1.x)*(p2.x-p1.x)+($y-p1.y)*(p2.y-p1.y))/sq(dist(p1.x,p1.y,p2.x,p2.y));
    float x = p1.x+u*(p2.x-p1.x);
    float y = p1.y+u*(p2.y-p1.y);
    if(x>p1.x&&x>p2.x){ x = max(p1.x,p2.x); }
    if(x<p1.x&&x<p2.x){ x = min(p1.x,p2.x); }
    if(y>p1.y&&y>p2.y){ y = max(p1.y,p2.y); }
    if(y<p1.y&&y<p2.y){ y = min(p1.y,p2.y); }
    Point p = new Point(x,y,0);
    return p;
  }
  public Point get_intersect_point(Segment $s){
    //xy coordinates only :: must test for on-segment intersection before using this, otherwise it will project the segments as lines to their intersection point :: it will also return {9999,9999,9999} for parallel lines (should be tested for prior to using this)
    float pbua1 = (p2.x-p1.x)*($s.p1.y-p1.y)-(p2.y-p1.y)*($s.p1.x-p1.x);
    float pbu2 = (p2.y-p1.y)*($s.p2.x-$s.p1.x)-(p2.x-p1.x)*($s.p2.y-$s.p1.y);
    float x = 9999;
    float y = 9999;
    if(pbu2!=0){
      float pbua = pbua1/pbu2;
      x = $s.p1.x+($s.p2.x-$s.p1.x)*pbua;
      y = $s.p1.y+($s.p2.y-$s.p1.y)*pbua;
    }
    Point p = new Point(x,y,0);
    return p;
  }
  public void echo(int $indent){
    String indent = "";
    while(indent.length()<$indent){
      indent += "  ";
    }
    println(indent+"---- SEGMENT #"+n+" ----");
    p1.echo($indent+1);
    p2.echo($indent+1);
  }
  public void render(){
    beginShape();
    vertex(p1.x,p1.y,p1.z);
    vertex(p2.x,p2.y,p2.z);
    endShape();
  }
}

//------ LINE ------//
class Line{
  int n;
  float a,b;
  Segment s;
  Line(Point $p1,Point $p2){
    n = global.nlines;
    global.lines[global.nlines++] = this;
    if($p2.x-$p1.x!=0){
      a = ($p2.y-$p1.y)/($p2.x-$p1.x);
    }else{
      a = 999999999;
    }
    b = $p1.y-a*$p1.x;
    s = new Segment($p1,$p2);
  }
  public void crop(){
    if(global.stage.is_splittable(this)){
      s = get_polygon_chord(global.stage);
    }
  }
  public boolean is_ray_intersecting(float $x,float $y,boolean $inclusive){
    //xy coordinates only :: for inside-outside testing, etc. :: ray points right, starting from the given coordinates
    float pbua1 = (s.p2.x-s.p1.x)*($y-s.p1.y)-(s.p2.y-s.p1.y)*($x-s.p1.x);
    float pbu2 = s.p2.y-s.p1.y;
    if(pbu2==0){ return false; }
    float pbua = pbua1/pbu2;
    if($inclusive){
      if(pbua>=0){ return true; }
    }else{
      if(pbua>0){ return true; }
    }
    return false;
  }
  public boolean is_line_intersecting(Line $l){
    //xy coordinates only
    float pbu2 = (s.p2.y-s.p1.y)*($l.s.p2.x-$l.s.p1.x)-(s.p2.x-s.p1.x)*($l.s.p2.y-$l.s.p1.y);
    if(pbu2==0){ return false; }
    return false;
  }
  public boolean is_segment_intersecting(Segment $s,boolean $inclusive){
    if($s.is_line_intersecting(this,$inclusive)){ return true; }
    return false;
  }
  public boolean is_coord_on(float $x,float $y){
    //xy coordinates only
    float d = get_distance_to($x,$y);
    if(d<=distance_tolerance){
      return true;
    }
    return false;
  }
  public Segment get_polygon_chord(Polygon $p){
    //xy coordinates only
    Point[] temp1 = new Point[$p.nsegments];
    int n1 = 0;
    for(int i=0;i<$p.nsegments;i++){
      if($p.segments[i].is_line_intersecting(this,true)){
        temp1[n1] = $p.segments[i].get_intersect_point(s);
        n1++;
      }
    }
    Point p1 = new Point(9999,9999,9999);
    Point p2 = new Point(9999,9999,9999);
    if(n1>=2){
      if(n1!=2){
        //remove duplicates in case of on-point intersection
        Point[] temp2 = new Point[n1];
        int n2 = 0;
        for(int i=0;i<n1;i++){
          boolean found = false;
          for(int j=0;j<n2;j++){
            if(temp1[i].get_distance_to(temp2[j].x,temp2[j].y,temp2[j].z)<distance_tolerance){
              found = true;
              break;
            }
          }
          if(!found){
            temp2[n2] = temp1[i];
            n2++;
          }
        }
        if(n2==2){
          p1.move_to(temp2[0].x,temp2[0].y,temp2[0].z);
          p2.move_to(temp2[1].x,temp2[1].y,temp2[1].z);
        }
      }else{
        p1.move_to(temp1[0].x,temp1[0].y,temp1[0].z);
        p2.move_to(temp1[1].x,temp1[1].y,temp1[1].z);
      }
    }
    Segment s1 = new Segment(p1,p2);
    return s1;
  }
  public Point get_closest_point(float $x,float $y){
    //xy coordinates only
    float u = (($x-s.p1.x)*(s.p2.x-s.p1.x)+($y-s.p1.y)*(s.p2.y-s.p1.y))/sq(dist(s.p1.x,s.p1.y,s.p2.x,s.p2.y));
    float x = s.p1.x+u*(s.p2.x-s.p1.x);
    float y = s.p1.y+u*(s.p2.y-s.p1.y);
    Point p = new Point(x,y,0);
    return p;
  }
  public Point get_intersect_point(Line $l){
    //xy coordinates only
    Point p = s.get_intersect_point($l.s);
    return p;
  }
  public Segment get_segment_to_line(Line $l){
    //must test for zero distance before using this function
    float p13x,p13y,p13z,p43x,p43y,p43z,p21x,p21y,p21z,d1343,d4321,d1321,d4343,d2121,pbmua,pbmub;
    p13x = s.p1.x-$l.s.p1.x; p13y = s.p1.y-$l.s.p1.y; p13z = s.p1.z-$l.s.p1.z;
    p43x = $l.s.p2.x-$l.s.p1.x; p43y = $l.s.p2.y-$l.s.p1.y; p43z = $l.s.p2.z-$l.s.p1.z;
    p21x = s.p2.x-s.p1.x; p21y = s.p2.y-s.p1.y; p21z = s.p2.z-s.p1.z;
    d1343 = p13x*p43x+p13y*p43y+p13z*p43z; d4321 = p43x*p21x+p43y*p21y+p43z*p21z; d1321 = p13x*p21x+p13y*p21y+p13z*p21z; d4343 = p43x*p43x+p43y*p43y+p43z*p43z; d2121 = p21x*p21x+p21y*p21y+p21z*p21z;
    pbmua = (d1343*d4321-d1321*d4343)/(d2121*d4343-d4321*d4321);
    pbmub = (d1343+d4321*pbmua)/d4343;
    Point p1 = new Point(s.p1.x+pbmua*p21x,s.p1.y+pbmua*p21y,s.p1.z+pbmua*p21z);
    Point p2 = new Point($l.s.p1.x+pbmub*p43x,$l.s.p1.y+pbmub*p43y,$l.s.p1.z+pbmub*p43z);
    return new Segment(p1,p2);
  }
  public float get_distance_to(float $x,float $y){
    //xy coordinates only
    float u = (($x-s.p1.x)*(s.p2.x-s.p1.x)+($y-s.p1.y)*(s.p2.y-s.p1.y))/sq(dist(s.p1.x,s.p1.y,s.p2.x,s.p2.y));
    float x = s.p1.x+u*(s.p2.x-s.p1.x);
    float y = s.p1.y+u*(s.p2.y-s.p1.y);
    return dist(x,y,$x,$y);
  }
  public void echo(int $indent){
    String indent = "";
    while(indent.length()<$indent){
      indent += "  ";
    }
    println(indent+"---- LINE #"+n+" ----");
    println(indent+"a: "+a+"  b: "+b);
    s.echo($indent+1);
  }
  public void render(){
    s.render();
  }
}

class Polyline{
  Segment[] segments = new Segment[1000];
  int nsegments = 0;
  int n;
  Polyline(Segment[] $segments){
    segments = $segments;
    nsegments = segments.length;
    n = global.npolylines;
    global.polylines[global.npolylines++] = this;
  }
  public float get_length(){
    float l = 0;
    for(int i=0;i<nsegments;i++){
      l += segments[i].get_length();
    }
    return l;
  }
  public void echo(int $indent){
    String indent = "";
    while(indent.length()<$indent){
      indent += "  ";
    }
    println(indent+"---- POLYLINE #"+n+" ----");
    for(int i=0;i<nsegments;i++){
      segments[i].echo($indent+1);
    }
  }
  public void render(){
    for(int i=0;i<nsegments;i++){
      segments[i].render();
    }
  }
}

//------ POLYGON ------//
class Polygon{
  int n,npoints,nsegments;
  Point[] points = new Point[1000];
  Segment[] segments = new Segment[1000];
  Polygon(Point[] $points){
    n = global.npolygons;
    global.polygons[global.npolygons++] = this;
    npoints = nsegments = 0;
    for(int i=0;i<$points.length;i++){
      add_point($points[i]);
    }
  }
  public void reset(){
    for(int i=0;i<nsegments;i++){
      segments[i].reset();
    }
  }
  public void add_point(Point $p){
    for(int i=0;i<npoints;i++){
      if($p.get_distance_to(points[i].x,points[i].y,points[i].z)<distance_tolerance){ return; }
    }
    points[npoints] = $p;
    npoints++;
    if(npoints>2){
      if(npoints==3){
        //create first closing segment
        segments[nsegments] = new Segment(points[npoints-2],points[npoints-1]);
        segments[nsegments+1] = new Segment(points[npoints-1],points[0]);
        nsegments += 2;
      }else{
        //replace closing segment
        segments[nsegments-1] = new Segment(points[npoints-2],points[npoints-1]);
        segments[nsegments] = new Segment(points[npoints-1],points[0]);
        nsegments++;
      }
    }else if(npoints>1){
      segments[nsegments] = new Segment(points[npoints-2],points[npoints-1]);
      nsegments++;
    }
  }
  public void delete_point(Point $p){
    //this is not an efficient way to remove points, just the easiest :: for reducing lots of large polygons, this should target only the segment to be removed, rather than rebuilding the entire polygon
    int ntemp = npoints;
    Point[] temp = new Point[1000];
    arraycopy(points,temp);
    points = new Point[1000];
    npoints = 0;
    nsegments = 0;
    for(int i=0;i<ntemp;i++){
      if(temp[i]!=$p){
        add_point(temp[i]);
      }
    }
  }
  public Polygon get_bounding_box(){
    //xy coordinates only
    float xmin = 9999; float ymin = 9999; float xmax = -9999; float ymax = -9999;
    for(int i=0;i<npoints;i++){
      xmin = min(points[i].x,xmin);
      ymin = min(points[i].y,ymin);
      xmax = max(points[i].x,xmax);
      ymax = max(points[i].y,ymax);
    }
    Point[] ps = new Point[4]; ps[0] = new Point(xmin,ymin,0); ps[1] = new Point(xmax,ymin,0); ps[2] = new Point(xmax,ymax,0); ps[3] = new Point(xmin,ymax,0);
    Polygon p = new Polygon(ps);
    return p;
  }
  public void simplify(){
    //xy coordinates only
    if(nsegments>2){
      for(int i=0;i<nsegments;i++){
        if(abs(segments[i].get_slope()-segments[(i+1)%nsegments].get_slope())<slope_tolerance){
          delete_point(points[(i+1)%npoints]);
          simplify();
          break;
        }
      }
    }
  }
  public boolean is_complex(){
    //xy coordinates only
    for(int i=0;i<nsegments;i++){
      for(int j=i+1;j<nsegments;j++){
        if(j==i+1||j==(i+nsegments-1)%nsegments){
          if(segments[i].is_segment_intersecting(segments[j],false)){ return true; }
        }else{
          if(segments[i].is_segment_intersecting(segments[j],true)){ return true; }
        }
      }
    }
    return false;
  }
  public boolean is_coord_inside(float $x,float $y,boolean $inclusive){
    //xy coordinates only
    int ints = 0;
    for(int i=0;i<nsegments;i++){
      if(segments[i].is_ray_intersecting($x,$y,$inclusive)){
        ints++;
      }
    }
    if(ints%2==0){ return false; }
    return true;
  }
  public boolean is_splittable(Line $l){
    //xy coordinates only
    if(is_complex()){ return false; }
    int ints = 0;
    for(int i=0;i<npoints;i++){
      if($l.is_coord_on(points[i].x,points[i].y)){
        ints++;
      }else if($l.is_segment_intersecting(segments[i],false)){
        ints++;
      }
    }
    if(ints==2){ return true; }
    return false;
  }
  public float get_distance_to(float $x,float $y){
    //xy coordinates only
    float d = segments[0].get_distance_to($x,$y);
    for(int i=1;i<nsegments;i++){
      float temp = segments[i].get_distance_to($x,$y);
      if(temp<d){ d = temp; }
    }
    return d;
  }
  public Point get_closest_point(float $x,float $y){
    //xy coordinates only
    float d = segments[0].get_distance_to($x,$y);
    Point p = segments[0].get_closest_point($x,$y);
    for(int i=1;i<nsegments;i++){
      float temp = segments[i].get_distance_to($x,$y);
      if(temp<d){
        d = temp;
        p = segments[i].get_closest_point($x,$y);
      }
    }
    return p;
  }
  public Point get_centroid_point(){
    float xs = 0; float ys = 0; float zs = 0;
    for(int i=0;i<npoints;i++){
      xs += points[i].x; ys += points[i].y; zs += points[i].z;
    }
    xs /= npoints; ys /= npoints; zs /= npoints;
    Point p = new Point(xs,ys,zs);
    return p;
  }
  public Polygon get_convex_hull(){
    //xy coordinates only :: returns the convex hull for this polygon's set of points :: this is a relatively costly function as it is scripted here
    int ip = 0; int ntemp = 0; Point[] temp = new Point[1000]; float[] as = new float[1000]; float[] ds = new float[1000];
    //find the topmost point
    for(int i=1;i<npoints;i++){
      if(points[i].y<points[ip].y){ ip = i; }
    }
    //if two points are colinear with the topmost point, take only the one furthest away
    for(int i=0;i<npoints;i++){
      if(i!=ip){
        float a = (atan2(points[ip].y-points[i].y,points[ip].x-points[i].x)+PI)%(TWO_PI);
        float d = dist(points[ip].x,points[ip].y,points[i].x,points[i].y);
        if(d>0){
          boolean found = false;
          for(int j=0;j<ntemp;j++){
            if(a==as[j]){
              if(d>ds[j]){
                //remove the closer point
                Point[] temp1 = (Point[]) subset(temp,0,j); Point[] temp2 = (Point[]) subset(temp,j+1,temp.length-(j+1)); temp = (Point[]) concat(temp1,temp2);
                float[] as1 = subset(as,0,j); float[] as2 = subset(as,j+1,temp.length-(j+1)); as = concat(as1,as2);
                float[] ds1 = subset(ds,0,j); float[] ds2 = subset(ds,j+1,temp.length-(j+1)); ds = concat(ds1,ds2);
                ntemp--; j--;
              }else{
                //skip this point
                found = true;
                break;
              }
            }
          }
          if(!found){
            //add point to array
            temp[ntemp] = points[i]; as[ntemp] = a; ds[ntemp] = d;
            ntemp++;
          }
        }
      }
    }
    //bubble sort the arrays according to angle
    for(int i=0;i<ntemp-1;i++){
      for(int j=0;j<ntemp-1-i;j++){
        if(as[j]>as[j+1]){
          Point temp1 = temp[j]; temp[j] = temp[j+1]; temp[j+1] = temp1;
          float as1 = as[j]; as[j] = as[j+1]; as[j+1] = as1;
        }
      }
    }
    temp[ntemp] = points[ip];
    ntemp++;
    temp = (Point[]) subset(temp,0,ntemp);
    temp = (Point[]) reverse(temp);
    //find convex hull points
    if(ntemp>3){
      for(int i=0;i<ntemp;i++){
        float pbside = (temp[(i+1+ntemp)%ntemp].y-temp[(i+ntemp)%ntemp].y)*(temp[(i+2+ntemp)%ntemp].x-temp[(i+ntemp)%ntemp].x)-(temp[(i+1+ntemp)%ntemp].x-temp[(i+ntemp)%ntemp].x)*(temp[(i+2+ntemp)%ntemp].y-temp[(i+ntemp)%ntemp].y);
        if(pbside<0){
          //remove the point if it lies on the wrong side
          if((i+1+ntemp)%ntemp==0){
            temp = (Point[]) subset(temp,1,ntemp-1);
          }else if((i+1+ntemp)%ntemp==ntemp-1){
            temp = (Point[]) subset(temp,0,ntemp-1);
          }else{
            Point[] temp1 = (Point[]) subset(temp,0,(i+1+ntemp)%ntemp);
            Point[] temp2 = (Point[]) subset(temp,(i+2+ntemp)%ntemp,ntemp-(i+2+ntemp)%ntemp);
            temp = (Point[]) concat(temp1,temp2);
          }
          i -= 2; ntemp--;
          temp = (Point[]) subset(temp,0,ntemp);
        }
      }
    }
    //return new polygon
    temp = (Point[]) subset(temp,0,ntemp);
    Polygon p = new Polygon(temp);
    p.simplify();
    return p;
  }
  public Polygon[] split(Line $l){
    //xy coordinates only :: must test to make sure it is a simple polygon first or you'll get strange results
    int ip = 0; Point[] ints = new Point[npoints]; float[] intsi = new float[npoints+2]; int nints = 0; Point[] temp1 = new Point[npoints+2]; float[] temp1i = new float[npoints+2]; int ntemp1 = 0; Point[] temp2 = new Point[npoints+2]; int ntemp2 = 0;
    for(int i=0;i<npoints;i++){
      if($l.is_coord_on(points[i].x,points[i].y)){
        ints[nints] = points[i]; intsi[nints] = i; nints++;
        temp1[ntemp1] = points[i]; temp1i[ntemp1] = i; ntemp1++;
        ip = i;
      }else if($l.is_segment_intersecting(segments[i],false)){
        Point lint = $l.s.get_intersect_point(segments[i]);
        ints[nints] = lint; intsi[nints] = i+0.5f; nints++;
        if(nints%2==0){ temp1[ntemp1] = points[i]; temp1i[ntemp1] = i; ntemp1++; }
        temp1[ntemp1] = lint; temp1i[ntemp1] = i+0.5f; ntemp1++; 
        ip = i;
      }else if(nints%2==1){
        temp1[ntemp1] = points[i]; temp1i[ntemp1] = i; ntemp1++;
      }
    }
    for(int i=ip;i<ip+npoints;i++){
      boolean found = false;
      for(int j=0;j<ntemp1;j++){
        if(i%npoints==temp1i[j]){ found = true; break; }
      }
      if(!found){ temp2[ntemp2] = points[i%npoints]; ntemp2++; }
    }
    temp2[ntemp2] = ints[0]; ntemp2++;
    temp2[ntemp2] = ints[1]; ntemp2++;
    temp1 = (Point[]) subset(temp1,0,ntemp1);
    temp2 = (Point[]) subset(temp2,0,ntemp2);
    Polygon p1 = new Polygon(temp1);
    Polygon p2 = new Polygon(temp2);
    Polygon[] ps = {p1,p2};
    return ps;
  }
  public void echo(int $indent){
    String indent = "";
    while(indent.length()<$indent){
      indent += "  ";
    }
    println(indent+"---- POLYGON #"+n+" ----");
    for(int i=0;i<nsegments;i++){
      segments[i].echo($indent+1);
    }
  }
  public void render(){
    beginShape();
    for(int i=0;i<npoints;i++){
      vertex(points[i].x,points[i].y,points[i].z);
    }
    endShape(CLOSE);
    /*
    for(int i=0;i<npoints;i++){
      points[i].render();
    }
    */
  }
}

class Plane{
  Vector nv;
  int n;
  float d;
  Point p1,p2,p3;
  Vector pv1,pv2,pv3,v12,v23,v31;
  Plane(Point $p1,Point $p2,Point $p3){
    p1 = $p1; p2 = $p2; p3 = $p3;
    n = global.nplanes;
    global.planes[global.nplanes++] = this;
    v12 = new Vector(p2.x-p1.x,p2.y-p1.y,p2.z-p1.z); v23 = new Vector(p3.x-p2.x,p3.y-p2.y,p3.z-p2.z); v31 = new Vector(p1.x-p3.x,p1.y-p3.y,p1.z-p3.z);
    nv = get_normal();
    d = get_d();
  }
  public Vector get_normal(){
    v12.vx = p2.x-p1.x; v12.vy = p2.y-p1.y; v12.vz = p2.z-p1.z; v23.vx = p3.x-p2.x; v23.vy = p3.y-p2.y; v23.vz = p3.z-p2.z;
    nv = cross_product(v12,v23);
    nv.unitize();
    nv.scale(100);
    return nv;
  }
  public float get_d(){
    d = -dot_product_1p(nv,p1);
    return d;
  }
  public int get_point_side(Point $p){
    float side = dot_product_1p(nv,$p)-d;
    if(side<0){ return 1; }
    if(side>0){ return -1; }
    return 0;
  }
  public boolean is_plane_parallel(Plane $pl1){
    Vector v = cross_product($pl1.nv,nv);
    if(v.m==0){ return true; }
    return false;
  }
  public boolean is_line_intersecting(Line $l){
    float pbmud = nv.vx*($l.s.p2.x-$l.s.p1.x)+nv.vy*($l.s.p2.y-$l.s.p1.y)+nv.vz*($l.s.p2.z-$l.s.p1.z);
    if(pbmud==0){ return false; } //plane and line do not intersect
    return true;
  }
  public boolean is_segment_intersecting(Segment $s,boolean $inclusive){
    float pbmud,pbmu;
    nv = get_normal();
    d = get_d();
    pbmud = nv.vx*($s.p2.x-$s.p1.x)+nv.vy*($s.p2.y-$s.p1.y)+nv.vz*($s.p2.z-$s.p1.z);
    if(pbmud==0){ return false; } //plane and line do not intersect
    pbmu = -(d+nv.vx*$s.p1.x+nv.vy*$s.p1.y+nv.vz*$s.p1.z)/pbmud;
    if($inclusive){ if(pbmu<0||pbmu>1){ return false; } }else{ if(pbmu<=0||pbmu>=1){ return false; } } //intersection is not on segment
    return true;
  }
  public Point get_3plane_intersection_point(Plane $pl1,Plane $pl2){
    //must check for parallel planes first
    Vector n12,n23,n31; float px,py,pz,pbden;
    nv = get_normal();
    d = get_d();
    n12 = cross_product($pl1.nv,$pl2.nv); n23 = cross_product($pl2.nv,nv); n31 = cross_product(nv,$pl1.nv);
    pbden = dot_product($pl1.nv,n23);
    px = ($pl1.d*n23.vx+$pl2.d*n31.vx+d*n12.vx)/pbden;
    py = ($pl1.d*n23.vy+$pl2.d*n31.vy+d*n12.vy)/pbden;
    pz = ($pl1.d*n23.vz+$pl2.d*n31.vz+d*n12.vz)/pbden;
    return new Point(px,py,pz);
  }
  public void echo(int $indent){
    String indent = "";
    while(indent.length()<$indent){
      indent += "  ";
    }
    println(indent+"---- PLANE #"+n+" ----");
    nv.echo($indent+1);
  }
}

class Face extends Polygon{
  Point p1,p2,p3;
  Plane pl;
  Face(Point[] $points){
    super($points);
    p1 = $points[0]; p2 = $points[1]; p3 = $points[2];
    pl = new Plane(p1,p2,p3);
  }
  public void orient(Point $p){
    if(pl.get_point_side($p)>0){
      flip();
    }
  }
  public void flip(){
    Point tp2 = p2;
    p2 = p3; p3 = tp2;
  }
  public boolean is_acute(){
    float d1,d2,d3,x,h,y,a1,a2,a3;
    d1 = p1.get_distance_to(p2.x,p2.y,p2.z); d2 = p2.get_distance_to(p3.x,p3.y,p3.z); d3 = p3.get_distance_to(p1.x,p1.y,p1.z);
    x = (sq(d1)+sq(d2)-sq(d3))/(2*d2); h = sqrt(sq(d1)-sq(x)); y = d2-x;
    a1 = (atan(h/x)+PI)%PI; a2 = (atan(h/y)+PI)%PI; a3 = PI-(a1+a2);
    if(a1<HALF_PI&&a2<HALF_PI&&a3<HALF_PI){ return true; }
    return false;
  }
  public Point get_incenter_point(){
    float d1,d2,d3; Point p;
    d1 = p2.get_distance_to(p3.x,p3.y,p3.z); d2 = p3.get_distance_to(p1.x,p1.y,p1.z); d3 = p1.get_distance_to(p2.x,p2.y,p2.z);
    p = new Point(((d1*p1.x)+(d2*p2.x)+(d3*p3.x))/(d1+d2+d3),((d1*p1.y)+(d2*p2.y)+(d3*p3.y))/(d1+d2+d3),((d1*p1.z)+(d2*p2.z)+(d3*p3.z))/(d1+d2+d3));
    return p;
  }
  public Point get_circumcenter_point(){
    Vector p12,p23,n12,n23; Point mid12,mid23; Line l12,l23; Segment s;
    pl.nv = pl.get_normal();
    p12 = new Vector(p2.x-p1.x,p2.y-p1.y,p2.z-p1.z); p23 = new Vector(p3.x-p2.x,p3.y-p2.y,p3.z-p2.z);
    n12 = cross_product(pl.nv,p12); n23 = cross_product(pl.nv,p23);
    mid12 = new Point((p2.x+p1.x)/2,(p2.y+p1.y)/2,(p2.z+p1.z)/2); mid23 = new Point((p3.x+p2.x)/2,(p3.y+p2.y)/2,(p3.z+p2.z)/2);
    l12 = new Line(mid12,new Point(mid12.x+n12.vx,mid12.y+n12.vy,mid12.z+n12.vz));
    l23 = new Line(mid23,new Point(mid23.x+n23.vx,mid23.y+n23.vy,mid23.z+n23.vz));
    s = l12.get_segment_to_line(l23);
    return s.p1;
  }
  public boolean is_line_intersecting(Line $l){
    float pbmud,pbmu,px,py,pz,pa1x,pa1y,pa1z,pa1m,pa2x,pa2y,pa2z,pa2m,pa3x,pa3y,pa3z,pa3m,a;
    pl.nv = pl.get_normal();
    pl.d = pl.get_d();
    pbmud = pl.nv.vx*($l.s.p2.x-$l.s.p1.x)+pl.nv.vy*($l.s.p2.y-$l.s.p1.y)+pl.nv.vz*($l.s.p2.z-$l.s.p1.z);
    if(pbmud==0){ return false; } //plane and line do not intersect
    pbmu = -(pl.d+pl.nv.vx*$l.s.p1.x+pl.nv.vy*$l.s.p1.y+pl.nv.vz*$l.s.p1.z)/pbmud;
    px = $l.s.p1.x+pbmu*($l.s.p2.x-$l.s.p1.x); py = $l.s.p1.y+pbmu*($l.s.p2.y-$l.s.p1.y); pz = $l.s.p1.z+pbmu*($l.s.p2.z-$l.s.p1.z);
    pa1x = p1.x-px; pa1y = p1.y-py; pa1z = p1.z-pz; pa1m = mag(pa1x,pa1y,pa1z); pa1x /= pa1m; pa1y /= pa1m; pa1z /= pa1m;
    pa2x = p2.x-px; pa2y = p2.y-py; pa2z = p2.z-pz; pa2m = mag(pa2x,pa2y,pa2z); pa2x /= pa2m; pa2y /= pa2m; pa2z /= pa2m;
    pa3x = p3.x-px; pa3y = p3.y-py; pa3z = p3.z-pz; pa3m = mag(pa3x,pa3y,pa3z); pa3x /= pa3m; pa3y /= pa3m; pa3z /= pa3m;
    a = acos(pa1x*pa2x+pa1y*pa2y+pa1z*pa2z)+acos(pa2x*pa3x+pa2y*pa3y+pa2z*pa3z)+acos(pa3x*pa1x+pa3y*pa1y+pa3z*pa1z);
    if(abs(a-TWO_PI)>angle_tolerance){ return false; } //intersection is not within face
    return true;
  }
  public boolean is_segment_intersecting(Segment $s,boolean $inclusive){
    float pbmud,pbmu,px,py,pz,pa1x,pa1y,pa1z,pa1m,pa2x,pa2y,pa2z,pa2m,pa3x,pa3y,pa3z,pa3m,a;
    pl.nv = pl.get_normal();
    pl.d = pl.get_d();
    pbmud = pl.nv.vx*($s.p2.x-$s.p1.x)+pl.nv.vy*($s.p2.y-$s.p1.y)+pl.nv.vz*($s.p2.z-$s.p1.z);
    if(pbmud==0){ return false; } //plane and line do not intersect
    pbmu = -(pl.d+pl.nv.vx*$s.p1.x+pl.nv.vy*$s.p1.y+pl.nv.vz*$s.p1.z)/pbmud;
    if($inclusive){ if(pbmu<0||pbmu>1){ return false; } }else{ if(pbmu<=0||pbmu>=1){ return false; } } //intersection is not on segment
    px = $s.p1.x+pbmu*($s.p2.x-$s.p1.x); py = $s.p1.y+pbmu*($s.p2.y-$s.p1.y); pz = $s.p1.z+pbmu*($s.p2.z-$s.p1.z);
    pa1x = p1.x-px; pa1y = p1.y-py; pa1z = p1.z-pz; pa1m = mag(pa1x,pa1y,pa1z); pa1x /= pa1m; pa1y /= pa1m; pa1z /= pa1m;
    pa2x = p2.x-px; pa2y = p2.y-py; pa2z = p2.z-pz; pa2m = mag(pa2x,pa2y,pa2z); pa2x /= pa2m; pa2y /= pa2m; pa2z /= pa2m;
    pa3x = p3.x-px; pa3y = p3.y-py; pa3z = p3.z-pz; pa3m = mag(pa3x,pa3y,pa3z); pa3x /= pa3m; pa3y /= pa3m; pa3z /= pa3m;
    a = acos(pa1x*pa2x+pa1y*pa2y+pa1z*pa2z)+acos(pa2x*pa3x+pa2y*pa3y+pa2z*pa3z)+acos(pa3x*pa1x+pa3y*pa1y+pa3z*pa1z);
    if(abs(a-TWO_PI)>angle_tolerance){ return false; } //intersection is not within face
    return true;
  }
  public boolean cull(){
    float p1x,p1y,p2x,p2y,p3x,p3y;
    p1x = screenX(p1.x,p1.y,p1.z); p1y = screenY(p1.x,p1.y,p1.z);
    p2x = screenX(p2.x,p2.y,p2.z); p2y = screenY(p2.x,p2.y,p2.z);
    p3x = screenX(p3.x,p3.y,p3.z); p3y = screenY(p3.x,p3.y,p3.z);
    if(((p2y-p1y)/(p2x-p1x)-(p3y-p1y)/(p3x-p1x)<0)^(p1x<=p2x==p1x>p3x)){ return false; }
    return true;
  }
  public void render(boolean $cull){
    if(!$cull||!cull()){
      super.render();
    }
  }
}

  static public void main(String args[]) {     PApplet.main(new String[] { "electronsphere" });  }}